Path Integral Approach to non-Markovian First-Passage Time Problems 
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The computation of the probability of the first-passage time through a given threshold of a 
stochastic process is a classic problem that appears in many branches of physics. When the stochastic 
dynamics is markovian, the probability admits elegant analytic solutions derived from the Fokker- 
Planck equation with an absorbing boundary condition while, when the underlying dynamics is non- 
markovian, the equation for the probability becomes non-local due to the appearance of memory 
terms, and the problem becomes much harder to solve. We show that the computation of the 
probability distribution and of the first-passage time for non-Markovian processes can be mapped 
into the evaluation of a path-integral with boundaries, and we develop a technique for evaluating 
perturbatively this path integral, order by order in the non-Markovian terms. 

PACS numbers: 05.40.-a, 02.50.Ey 



The computation of the statistical distribution of the 
times at which a stochastic process £(t) first reaches a 
given threshold (the so-called first-passage time problem) 
is a classic problem that appears in many different con- 
texts in physics, chemistry and biology. It is relevant for 
instance to problems appearing in reaction rate theory, 
micleation theory or neuron firing, to name just a few, 
and it is treated in a number of textbooks [l|, 0, 0] and 
reviews Q. When the underlying dynamics is marko- 
vian, the function Ii{x,t) that gives the probability dis- 
tribution that had a value x, in the continuum limit 
satisfies a Fokker-Planck (FP) equation. The fact that 
one is interested in the first-passage time problem means 
that one wants to discard the trajectories once they have 
reached for the first time a given threshold x c . This is 
implemented imposing on the FP equation an absorbing 
barrier boundary condition Tl(x c ,t) = 0. The FP equa- 
tion with this boundary condition, together with the ini- 
tial condition H(x,t = 0) = 6r>(x = 0), where 80 is the 
Dirac delta, can be elegantly solved using the method of 
images, and one finds [fj 



IL(x;t) 



1 



e -*7(2t) _ e -(2x c -x) 2 /(2t) 



(1) 



When the underlying dynamics is non-Markovian, how- 
ever, the problem becomes much more difficult. The 
system acquires memory properties and the probability 
IL(x, t) no longer satisfy a simple diffusion equation such 
as the FP equation. Furthermore, the correctness of the 
"absorbing barrier" boundary condition is now far from 
obvious @,0|- For these reasons, first-passage problems 
for non-Markovian processes are known to be very hard 
to solve, and have been attacked in various way, see e.g. 
d, E3, EH E3, Ell an d references therein. Our orig- 
inal interest in the problem arose from a specific ques- 
tion in cosmology, namely the computation of the mass 
distribution of dark matter halos generated by the evo- 
lution of non-Gaussian primordial density fluctuations, 



which can indeed be formulated as a first-passage time 
problem with non-Markovian dynamics [14| . We think 
however that the techniques that we have developed in 
refs. [H, El, E3 , and which allowed us to solve our prob- 
lem, have a broader interest, and we find it useful to 
present them here in a more general context. 

Let £(t) be a variable that evolves stochastically with 
time t, with (£(t)) = 0. We consider an ensemble of 
trajectories starting at to = from an initial position 
£(0) = xo, and we follow them for a time t. We discretize 
the interval [0, t] in steps At — e, so tk — ke with k = 
X,..,n. A trajectory is then defined by the collection 
of values {xi, . . . ,x n }, such that = Xk- There is 

no absorbing barrier, i.e. £(t) is allowed to range freely 
from —00 to +00. The probability density in the space 
of trajectories is 



W(x ;xi, . . .,x n ;t n ) = (JI<fo(£(*i) - Xi)) . (2) 
In terms of W we define 



i=i 



U e (x ;x n -t n ) = Yl / dxiW(x ; X\ 3 . . . 3 X n \ tjij , (3) 

i=l "'- 00 

where t n — ne = t and we will often write x n simply as 
x. So, M t (xo\x;t) is the probability density of arriving 
in x at time t, starting from xo at time to = 0, through 
trajectories that never exceeded x c . Observe that the 
final point x ranges over —00 < x < 00. For later use, we 
find useful to write explicitly that II depends also on the 
temporal discretization step e. We are finally interested 
in its continuum limit, II e= o, an d we will see in due course 
that taking the limit e — ► of LI e is non-trivial. 

The usefulness of II e is that it allows us to compute 
the first-crossing rate from first principles, without the 
need of postulating the existence of an absorbing barrier. 
Simply, the quantity J_° dxU e (xo;x;t) gives the prob- 
ability that at time t a trajectory always stayed in the 



2 



region x < x c , for all times smaller than t. The rate of 
change of this quantity is therefore equal to minus the 
rate at which trajectories cross for the first time the bar- 
rier, so the first-crossing rate is 



T{t) 



dx n d t lJe(x ;x n ;t) . 



(4) 



Observe that no reference to a hypothetical "absorbing 
barrier" is made in this formalism. We will see below 
how an effective absorbing barrier emerges from this mi- 
croscopic approach. 

The probability density W can be expressed in terms 
of the connected correlators . . . £i p ) c as [l[ 



W(x ; xx, . . . , x n ; t n ) = VX e 



2 L-ii=\ 



(5) 



n n 



p=2 



i\— 1 i p — 1 



where / VX = §± . . . ^ and £ = T he prob- 

lem is therefore reduced to computing the path-integral 
<j3j> , over variables Xj, bounded by x c , with W given by 
eq. ©. 

We first consider the simple case in which £ has gaus- 
sian statistics (so only the two-point connected func- 
tion is non- vanishing) , and obeys a Langevin equation 
£ = Tj(t) with a noise r/ whose correlator is a Dirac delta, 
(r/(t)r](t')) = S D (t — t'). In this case the 2-point correlator 
is easily computed 

(£(<0e(*i)>c = r dt T dt'{r){t) V (t')) = mhxitutj) , 
Jo Jo 

(6) 

and from this, performing the gaussian integrals in 
eq. (J5J), wc find 



W sau (x ;xx,...,x n ;t n ) 



1 



(27re)"/ 2 



(7) 

where we denote by VF gau the value of W when £ has 
gaussian statistics and obeys a Langevin equation with 
Dirac-delta noise. Not surprisingly, in this case we got 
the Wiener measure. To compute IIf au by performing 
directly the integrals over xx, ■ ■ ■ , x n -x in eq. ([3]), and 
then taking the limit e — * is very difficult, since the 
integrals in eq. ((3|) run only up to x c , and already the 
inner integral gives an error function whose argument 
involves the next integration variable. In [l5| we have 
then followed a different route. Using the explicit form 
we proved that 

/>OC 

nf u (x ; x; t+e) = / d(Ax) * e (Ax)nf "(x ; x-Ax; t), 

J x — x c 

(8) 

where * e (Aa;) = (27re)- 1/2 exp{-(Aa;) 2 /(2e)}. Equa- 
tion {5| generalizes the Chapman-Kolmogorov equation, 
to which it reduces if we send x c — > oo, i.e. if the 



integrations in eq. ((3|) are not bounded, and expresses 
the fact that the evolution corresponding to a Langevin 
equation with Dirac-delta noise is a markovian pro- 
cess. Equation ijHJ) allows us to compute the continuum 
limit of IIf au , as follows. In the limit e — > we have 
^ e (Ax) — > Sd(Ax). If x — x c < 0, the integral in eq. © 
includes the support of the Dirac delta, and we just get 
the trivial identity that nff^iro; x; t) is equal to itself. 
However, if x — x c > 0, the right-hand side vanishes 
and we get H^ Q (xo;x;t) — 0. Therefore we find that 
H^ (xo;x;t) = if x > x c . (If x = x c only one half 
of the support of ^ e is inside the integration region, so 
we get nf!lo(2:o; x c ; t) = (l/2)IIf*o(a;o;x c ; t), which again 
implies Tl^ (xo; x c ; t) = 0). This is the boundary con- 
dition that in the usual treatment is just imposed by 
hand, while here it follows from the formalism. Con- 
sider now eq. ||5J) when x < x c . In this case the zeroth- 
order term in e gives a trivial identity. Pursuing the 
expansion to higher orders one finds that, in the limit 
(x c — x)j\fi — * + , and therefore when x is fixed and 
strictly smaller than x c while e — > + , the dependence on 
the index e in LIf au can be expanded in integer powers of 



Ilf "(x ; x; t) = n£ u (a;o; x; t)+elJ^(x ;x; t) 



(9) 



Collecting terms of the same order in e we then find that, 
for x < x c , nflo(a;o; x;t) satisfies a FP equation. We 
therefore end up with a FP equation with the boundary 
condition Il^ (xo',x — x c ;t) — 0, so we recover eq. ([lj. 
We have therefore succeeded in deriving this standard 
result from our path integral approach. Observe that the 
boundary condition H^ q (xq; x = x c ; t) — emerges only 
when we take the continuum limit, and does not hold for 
finite e. 

Having computed the path-integral in the markovian 
case, we can tackle the problem of non-Markovian dy- 
namics treating the non-Markovian terms as perturba- 
tions. For illustration, we discuss the case of colored 
gaussian noise, i.e. again the only non-vanishing con- 
nected correlator is the two-point correlator, but now 
we take it to have the form (£(fi)£(fj)) c = min (tj,tj) + 
A(ti,tj), for some function A(ti,tj) = Ay. So we want 
to compute 



,(10) 



Tl e (x n ;t n ) = / dxx...dx n -x I V\ 



x exp < iX 



-[min(U,tj) + A(t t , tj)]XiXj 



where, for simplicity, wc set xq = and we eliminated 
it from the list of variables on which n e depends. We 
assume that A^ is proportional to a small parameter, 
and we expand perturbatively in Ay. Using Xke lXx = 
—id x e lXx and writing d/dxi = di, the first-order correc- 
tion to n e , that we denote by n^ 1 , is 



n e (x n , t n ) 



dxx.-.dxn-xdidjW^ . 



(11) 
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We rewrite the term Aijdidj separating explicitly the 
derivative d n = d/dx n from the derivatives di with i < n. 
Let us at first consider the case in which A(tj, tj) vanishes 
at least as tj — U as t{ — > tj. This was indeed the case 
in the application to the cosmological problem discussed 
in [l5j |. The more general case will be discussed later. In 
this case, using also = Aji, 



Y n Tl — l 



(12) 



i< 3 



En— 2 ^^n— 1 
i— 1 ^-^7 



When inserted into 



where Y^<j - z^=i ^j=i+v 

eq. (fTTj) the term X^fcTi 1 ^indid n brings a factor ^ i that, 
in the continuum limit, produces an integral over an in- 
termediate time U. Because of this dependence on the 
past history, we call this the "memory term" . Similarly, 
the double sum in X^i<j ^-ijdidj gives, in the continuum 
limit, a double integral over intermediate times ti and 
tj, and we call it the "memory-of- memory" term. Thus, 
nf 1 = nf cm + nf cm - mcm , where 



IT 



i=l 



J2^nd n dx 1 ...dx n -id i W e * a , (13) 



n 



mem — mem 



i<j 



= ^2 Ay / dxi . . . dxn-ididj W g 



(14) 

To compute the memory term we integrate di by parts 
and we make use of the fact that W gau satisfies 

I^ gau (x ;xi, . ..,Xi = x c , ... ,x n ;t n ) (15) 
= W sau (x ;x 1 ,...,x i = x c ;t i ) 

as can be checked from the explicit expression ([?]) ■ so 



dx 1 ...dx„_ 1 a i W gau 

= IIf» u (a;o;a: c ;t i )IIf u (a; c ;i n ;t n -ti). (16) 

In the continuum limit (if the integral converges, as we 
will check in a moment) , we replace Y^=i by (1/e) / " dU 
and we get 



n e= o {x n ',t n ) — d n I dti A(tj,t„) 

Jo 

x lim ~nf u {x )X c -,ti)TS.f u {x c -,x n -,t n -U) 
c— »o+ e 



(17) 



It is quite interesting to observe that the memory term is 
determined by the finite-e corrections to the markovian 
term. We found above that H^ q {xq\ x n ;t) vanishes for 
x n = x c . However, we see from eq. (|17l) that it is not 
enough to know that IIf au (xo; x c ; t) = for e — ► 0, but 
we also need to know how fast it goes to zero with e. 



For x — x c fixed and strictly negative, in the limit e — > 
we have seen above that the correction to n ga Q are 0(e). 
However, in eq. (fT7|) we need nf au for x = x c . In this case 
the form of the correction changes qualitatively. Techni- 
cally this comes from the fact that, after changing the in- 
tegration variables from Xi to = Xi/y2e, which makes 
the exponential factors in eq. {7J independent of e, the 
lower integration limit in eq. ([5]) becomes (x — x c )/\/2e. 
For x fixed and strictly smaller that x c , when e — > + 
this lower limit goes to — oo, while for x — x c it is zero 
for all e, resulting in a different form of the solution. The 
passage between the two regimes takes place when the 
lower limit of the integral is 0(1), i.e. when x — x c ~ \fe. 
So, in the continuum limit, the path integral in eq. (J3]) 
has three different regimes: for x strictly smaller than 
x c , with x c — x finite, it approaches eq. ([TJ), phis cor- 
rections 0(e). For x > x c and x — x c finite it is equal 
to zero, plus corrections which can be shown to be ex- 
ponentially small in e, (exp{— (x c — x) 2 /(2e)}). These 
two regimes are connected by an infinitesimal boundary 
layer, of thickness \x — x c \ — 0(y/e), where the correc- 
tions to the continuum solution are themselves 0(y/e). In 
particular, when x — x c , we find in [lil ] that, for xq < x c , 



n^ u (x ;x c ;t) = yfe 



^-x«) 2 /(2t) +0 (e). (18) 



We see that the factor y/e from nf au (xo; x c ; U) and a 
similar factor y/e from IIf au (a; c ; x n ; t n — ti), cancel the 
factor 1/e which comes from transforming the sum into 
an integral. The limit in eq. (jTTJ) is then finite, and 

[ tn dUA(U,t n ) x f x °~ x J 



rrmcm / . . \ _ 
ll e— K^n^nJ 

7T 



x exp 



~2U 



(x c 



(19) 



Recall that IIf au (a;o; x c ; ti) represents the probability 
density for trajectories that start at xo at the initial 
time and arrive at x c at time ti, staying always in the 
region x < x c , for a variable obeying gaussian statistics 
and driven by a Dirac delta noise. Thus, eq. (jTTJ) has a 
vivid diagrammatic interpretation in terms of a sum over 
the markovian trajectories that start at xo, touch for the 
first time the boundary x c at an intermediate time ti, but 
rather than crossing the threshold go back into the region 
x < x c , finally reaching a value x n at time t n . We see 
that the probability density IIf au (a;o; x c ; U) plays the role 
that in the perturbative expansion of the path integral 
in quantum field theory is played by the free propagator, 
and the whole complexity of non-Markovian dynamics 
enters through the presence of the boundary at x — x c . 

The memory-of-memory term can be computed in the 
same way. Now we integrate by parts the two derivatives 
didj in eq. (|14p . This leaves us with W evaluated in 
Xi = x c and Xj = x c . Using eq. (|15p we write it as a 
product of three terms, and 



II 



mem — mem / 



t n ) = ^ A y nf au ( 



i< 3 



4 




FIG. 1: The first crossing rate J-{t) from eq. ((22} for the non- 
Markovian case with k = 0.44 (blue solid line) compared to 
the markovian result obtained setting k — 0. In both cases 
we set x c = 1. 



x 11^ (xc^Xctj £2)11^ (x c ,x n ,t n tj] . 



(20) 



Diagrammatically, this corresponds to a sum over the 
markovian trajectories that start at xq, touch for the 
first time the boundary x c at an intermediate time ti, go 
back into the region x < x c , touch again the boundary at 
time tj, and finally reach the value x n at time t n , always 
staying in the region Xi < x c for i < n. The explicit 
computation shows that Hf au (x c ; x c ;t) — e/(v / 2~7r i 3 / 2 ). 
Then again the continuum limit is finite, and 



n™- mcm (^ = 0;x„;(„) 



1 



x / dti / dtj^-j- 

x cxp <^ - — - ______ j { ) 

Equations (flT)|) and (f_T|) provides an analytic expres- 
sion for the first-order non-Markovian corrections due 
to a non-trivial two-point function, under the assump- 
tion that A(ti,tj) vanishes, at least linearly, as U — > tj. 
Given the probability distribution, the first-crossing rate 
is given by eq. ((4]), and from this one can derive all the 
properties concerning the statistics of the time at which 
the threshold is first crossed. 



As an application, we illustrate our results setting 
Aij = Kti(tj — ti)/tj for ti < tj. (The value for ti > tj 
is obtained by symmetry, Ay- = Aji). This was indeed 
the form of Ajj in the problem studied in [l5l ]. where 
k ~ 0.44 played the role of the expansion parameter. All 
integrals can then be computed analytically, and for the 
first-crossing rate we find 



2^ t 3 / 2 



c/(2t) 



+ 



x 2 
2t 



(22) 

where T(0, z) is the incomplete Gamma function. The 
result is shown in Fig. [TJ together with the markovian 
case, which is obtained setting k = 0. 

In the above example, A(ti,tj) goes to zero linearly 
as tj — tj — > 0. If however A(ti,tj) goes to a non-zero 
constant at ti — tj, we see that the integral over dtj in 
eq. (|2Tj) diverges at the lower limit tj = ti. At the same 



time, in eq. (fl~2"l) we must also include a term propor- 
tional to An&f, since An is non-zero, and this term also 
leads to a divergent integral. After regularizing the sums 
over i,j one can extract the divergent part of both in- 
tegrals, which are both proportional to plus the 
finite part. The divergent parts of these two terms can- 
cel among them, and we remain with a finite result. We 
refer the reader to appendix B of ref. [l5[ for details. 

The same strategy can be applied to all higher-order 
correlators. In ref. [I?) we performed the computation 
expanding eqs. and © to linear order in the three- 
point correlator (Ci£,j£,k)- The result fits quite well the 
outcome of cosmological iV-body simulations with non- 
Gaussian initial conditions [l8| , giving further confidence 
in our technique. 

In conclusion, we have developed a very general 
method for computing systematically the non-Markovian 
contributions to the first-crossing rate, whenever they 
can be treated as perturbations of the markovian dy- 
namics. 
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